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1.  INTRODUCTION 


Composite  material  offers  a  great  potential  and  flexibility  in  structural  design  because  of  the  anisotropy 
of  material  properties,  unique  ply-by-ply  constructions,  and  novel  fabrication  methods.  To  extract  the 
maximum  performance,  structures,  in  general,  are  designed  with  various  ply  orientations  and  stacking 
sequence  from  layer  to  layer  within  the  structures.  This  flexibility  does  enhance  the  structural  design; 
however,  it  also  increases  the  degree  of  difficulty  for  analysis,  particularly  by  using  the  finite  element 
method  (FEM).  This  is  especially  true  for  a  thick-section  composite  structure,  which  may  consist  of 
thousands  of  anisotropic  layers. 

There  are  two  analytical  approaches  that  can  be  used  for  analysis  of  composite  structures  by  FEM, 
a  layer-by-layer  analysis  or  a  property  smearing  model.  The  layer-by-layer  analysis  approach  will,  in 
general,  result  in  a  huge  finite  element  model  with  thousands  of  elements  required  to  maintain  a  proper 
aspect  ratio  of  elements.  Therefore,  a  tremendous  computational  effort  is  required,  especially  for  a 
dynamic  analysis.  For  a  thick-section  large  composite  structure,  the  layer-by-layer  approach  is  not 
practical. 

Another  approach  is  to  use  the  smeared  (effective)  properties  for  the  elements.  Accordingly,  each 
element  consists  of  several  layers  and  material  blocks.  The  properties  of  elements  are  calculated  from  the 
properties  of  the  contained  layers  based  on  certain  assumptions.  The  effective  properties  of  the  input 
model  are  crucial  to  the  accuracy  of  FEM  analysis.  Several  models  based  on  either  "laminated  plate 
theory"  or  "rule  of  mixture"  have  been  developed  to  compute  the  effective  properties  for  use  for  FEM 
analysis.  However,  these  approaches  cannot  be  used  to  calculate  the  effective  properties  accurately  for 
a  very  common  element  with  irregular  geometry.  For  example,  a  taper-shaped  element  from  a  filament- 
wound  cone  model  is  illustrated  in  Figure  1. 

Enie  and  Rizzo  (1970),  Pagano  (1974),  and  Christensen  and  Zywicz  (1989)  derived  the  effective 
properties  from  laminated  plate  theory.  Particularly,  Pagano’s  model  is  an  exact  three-dimensional  (3-D) 
solution  calculated  from  a  laminated  plate.  These  models  were  all  developed  by  assuming  a  finite 
thickness  in  the  transverse  direction  (2-D  geometry).  In  general,  a  constant  interlaminar  shear  stress 
distribution  through  the  thickness  is  assumed  for  these  models.  Properties  calculated  from  these  models 
may  be  suitable  for  thin-shell  structures,  but  are  not  proper  for  a  thick  structure  or  a  block  element  since 
these  models  do  not  correctly  account  for  the  properties  in  the  transverse  direction,  especially  for  the 
transverse  shear  and  shear  coupling  properties.  In  addition,  the  plate  theory  is  limited  to  a  rectangular 
geometry  and  can  never  be  used  for  an  arbitrarily  shaped  geometry  which  commonly  occurs  in  FEM 
modeling. 
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Finite  element  model 


Figure  1.  Taper-shaped  elements  in  an  FEM  cone  model. 

Models  developed  by  Chou,  Carleone,  and  Hsu  (1971)  and  Sun  and  Li  (1988)  assumed  a  uniform 
displacement  in  the  planes  parallel  to  the  laminate  and  a  uniform  stress  in  the  transverse  direction  of  the 
plane.  The  transverse  (interlaminar)  shear  components  are  calculated  on  the  basis  of  volume  average.  The 
in-plane  properties  resulting  from  this  approach  are  correct.  However,  the  transverse  shear  properties  are 
not  accurate  enough  and,  thus,  not  suitable,  for  thick-section  structures,  which  generally  have  large 
transverse  shear  deformations. 

One  of  the  common  shortcomings  of  these  models  is  the  limitation  of  geometry.  The  laminate  is 
generally  restricted  to  uniform  thickness  flat  plate  or  thin  shell  configurations  with  layers  aligned  along 
element  boundaries.  Element  faces  must  be  rectangular  in  both  the  plane  of  the  laminate  and  the  through¬ 
thickness  directions.  This  limitation  makes  it  very  difficult  to  model  regions  containing  ply  drop-offs,  or 
layer  terminations,  and  impose  restrictions  on  the  capability  to  generate  finite  element  meshes  for  complex 
geometries  that  require  changes  in  mesh  density  and/or  arbitrarily  shaped  elements  that  cannot  be  readily 
aligned  with  the  layers  of  the  laminate. 

Figure  2  shows  a  generalized  case  of  an  element  in  a  region  containing  several  materials.  The 
effective  properties  of  the  element  certainly  cannot  be  calculated  correctly  with  any  of  the  models 
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Laminated  Block 


Figure  2.  Arbitrarily  shaped  element  with  multiple  anisotropic  layers  and  materials. 


mentioned  previously.  Commercial  packages  typically  use  the  volume  average  approach  for  computing 
effective  properties  in  elements  such  as  this,  leading  to  potential  inaccuracies  in  results,  especially  for 
irregular  shaped  elements.  Accordingly,  there  is  a  strong  need  to  develop  an  accurate  property  model  for 
FEM  applications. 


The  objective  of  this  investigation  is  to  develop  a  model  that  provides  accurate  3-D  effective  properties 
for  arbitrarily  shaped  solid  continuum  elements  containing  multiple  layers  and  materials  with  various 
orientations  and  shapes.  The  second  objective  will  be  to  develop  a  pre-processor  incorporating  the 
effective  property  model  in  generating  accurate  finite  element  representations  of  3-D  laminated  material 
structures  for  ABAQUS  and  DYNA3D. 


2.  MODEL  DEVELOPMENT 


Consider  a  portion  of  a  material  block  contained  within  some  internal  region  of  an  element’s  volume. 
Relative  to  the  global  frame  with  coordinates  (x1,  x2,  x3)  illustrated  in  Figure  3,  the  generalized  elastic 
constitutive  behavior  of  each  material  block  is  described  by  the  relations 


Figure  3.  Coordinates  transforms  of  material  blocks  and  elements. 

i  k  1 1 

where  tj  is  the  stress  tensor,  Ej  is  the  strain  tensor,  and  Cjk  represents  the  material  stiffness  tensor 

relative  to  the  global  frame. 

We  assume  the  deformation  field  within  the  volume  portion  contained  in  the  element  bounds  is 
continuous  such  that  it  can  be  approximated  as  a  function  of  the  displacements  at  the  comer  points  that 
bound  the  volume  portion. 

We  first  incorporate  the  transformation  of  global  coordinates  to  isoparametric  coordinates  within  the 
volume  portion 

xi=NaXai,  (1) 

where  the  Xj  are  the  global  coordinates  of  the  comer  points. 
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The  notation  used  herein  follows  the  standard  summation  convention  for  repeated  indices.  English 
indices  correspond  to  the  range  (1,2,3),  referring  to  the  three  independent  spacial  components  of  a  variable 
while  Greek  indices  correspond  to  the  (1,2,.. .T)  discrete  point  values  and  N^g1,  g2,  g3)  isoparametric 
interpolation  functions  used  for  describing  the  continuous  variation  of  a  parameter  within  a  volume.  The 
interpolation  functions,  written  in  terms  of  the  isoparametric  coordinates  g*  of  the  volume,  are  expressed 
as: 


Linear  formulation  (T=8): 


N^id-g1)  (1-g2)  (1-g3) 

O 

N2=i(l+g1)  (1-g2)  (1-g3) 

O 

N3=A(l+g1)  (l+g2)  (1-g3) 

O 

N4=I(l-g1)  (l+g2)  (1-g3) 

O 


N^J-d-g1)  (1-g2)  (l+g3) 

O 

N6=-l(l+g1)  (1-g2)  (l+g3) 

O 

N7=I(l+g1)  (l+g2)  (l+g3) 

O 

N^^l-g1)  (l+g2)  (l+g3) 


Quadratic  formulation  (T=20): 


N^-Id-g1)  (1-g2)  (1-g3)  (2+g1+g2+g3) 

O 

N2  =  -l(l+g1)  (1-g2)  (1-g3)  (2-g1+g2+g3) 

O 

N3  =  -i.(l+g1)  (l+g2)  (1-g3)  (2-g1-g2+g3) 

O 

N4=-i.d-gJ)  (l+g2)  (1-g3)  (2+g1-g2+g3) 

O 

N5— Id-g1)  (1-g2)  (l+g3)  (2+g  1+g2-g3) 

O 

N6=-i.(l+g1)  (1-g2)  (l+g3)  (2-g1+g2-g3) 

O 

N7— Id+g1)  (l+g2)  (l+g3)  (2-g1-g2-g3) 

O 

N^-id-g1)  (l+g2)  (l+g3)  (2+g1-g2-g3) 

O 
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N9=i-(l-g1)  (Ug1)  (1-g2)  (1-g3) 
4 

N10=l(l-g2)  (Ug2)  (Ug1)  (1-g3) 
4 

N^la-g1)  (Ug1)  (1+g2)  (1-g3) 
4 

N12=i.(l-g2)  (Ug2)  (l-g1)  (1-g3) 
4 

N^-id-g1)  (Ug1)  (1-g2)  (1+g3) 
4 

N14=i.(i-g2)  d+g2)  (1+g1)  d+g3) 

4 


N15=i(l-g1)  (Ug1)  (Ug2)  (Ug3) 
4 

N16=I(l-g2)  (l+g2)  (1-g1)  (l+g3) 
4 

N17=i-(l-g3)  (l+g3)  (1-g1)  (1-g2) 
4 

N18=i.(l-g3)  (l+g3)  (l+g1)  d-g2) 
4 

N19=i(l-g3)  (Ug3)  (l+g1)  (l+g2) 
4 


N20=i.(l-g3)  (Ug3)  (1-g1)  (Ug2) 
4 

Displacements  within  the  volume  portion  are  given  by 


u1  =  N°  U0\ 


(2) 


where  the  uj  represents  displacement  components  of  the  comer  points.  The  strain-displacement  relations 
are  thus 


•5-4 


S' 

+ 

^  ro 

i  aNa 

'agm  ui+  dgm 

[dx  j  dx  1 J 

2  dT” 

axj  a  ax 1 

'a 


(3) 


The  strain  energy  contributed  by  the  material  block  is  thus  expressed  as 


e(M)  =  - 


i 

dl 

r  dNa  3nP  dgm  dgn 

2 

ik 

(M) 

^  dg  m  3g  n  3x  1  3x 1 

(M) 


uj  Upk, 


(4) 


the  integration  being  performed  over  the  volume  portion  of  the  material  block  (M)  contained  within  the 
element  This  defines  the  stiffness  of  the  material  block  portion  as 


9Na  0NP  dg m 
dg  m  dg  n  dx  1 


dg" 

dx1 


dv 


l(M) 


(5) 
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so  that  the  strain  energy  can  be  written  as 


(M>  =  "2 


K 


“P 


ik 


uJ  uBk. 

(M)  a  P 


(6) 


The  components  of  force  contributed  by  the  material  block  at  comer  point  a  are  obtained  from 


®<M> 


K 


ap 

ik 


J(M) 


(7) 


Summing  over  all  material  blocks  that  are  common  to  comer  point  a,  results  in  an  expression  for  the  net 
external  applied  force  on  that  point 


Fi“  =  £ 

M 


Lk*J<m)  ' 


aP 


UR 


(8) 


Assuming  a  total  of  T  comer  points  contained  within  the  volume  of  the  element,  there  will  be  3  x  r 
equations  of  the  form  (8).  Let  the  first  £2  of  these  comer  points  correspond  to  the  nodes  of  the  element 
(£2  <  T),  the  next  A  correspond  to  comer  points  lying  on  the  surface  of  the  element  at  locations  other  than 
the  element  nodes  (0  <  A  <  (T-£2)),  and  the  last  (T-£2-A)  correspond  to  points  falling  within  the  interior 
of  the  element  boundary. 

If  we  represent  the  material  stiffness  of  the  overall  element  as  an  equivalent  homogeneous  anisotropic 
material  with  stiffness  tensor  C^,  then  the  total  strain  energy  for  the  element  is  given  by 


E  =  .L  C? 


ik 


ras* 

asp 

3q  m 

al”dv" 

dq  m 

dq  n 

ax  j 

ax1 

JV 


(9) 
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where  y,p=l,...Q  and  the  integration  is  over  the  entire  volume  of  the  element,  V;  SY  represents  the  element 
deformation  shape  functions;  and  qm  are  the  isoparametric  coordinates  of  the  transformation 

x^S^x/  SY  =  Sy(ql,  q2,  q3),  (10) 

with  Xy1  corresponding  to  the  global  coordinates  of  the  node  point  y.  The  force  applied  to  node  y  is  given 
by 


9SY  3SP  3q_^  dy 
3q  m  8q  n  Bx  i  8x  1 


(11) 


where  y,p=l,...,£2. 

To  establish  the  equivalent  material  stiffness,  we  equate  the  expression  for  total  strain  energy  in  the 
element  (9)  to  the  sum  of  the  strain  energies  contributed  by  the  individual  material  blocks  (6),  i.e.. 


Cjl 

'■'ik 


rasY 

asp 

aqm 

^”dv 

U  *  Upk  = 

/aqm 

3qn 

dx * 

dx  1 

Y  P 

(12) 


E  {  IKfldi)  UT  VP  *  2  IKi]f](M)  U,‘  Uk  ♦  [K^  Up‘  U„k  }. 

M 


where  y,p=l,...,Q  and  p,ro=(£2+l),...,r.  For  convenience,  we  define  the  following: 


(  «T  asp  dq m  aq° 

^  ^  dq  m  3q  n  3x  1  3x  1 


BiP  -  E 

M 


K;' 


HP 


ik 


(M) 


(13a) 


(13b) 
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Using  this  notation,  we  can  write  (12)  as 


c i  A»  U,‘  Upk  =  Bj  u;  Upk  *  2  B*  UT‘  Up1  *  b|“  Up*  U„k, 


(14) 


where  y,p=l . Q  and  p,<n=(Q+l),...,r. 

As  will  be  shown,  the  displacements  of  comer  points  p=(£2+l),...,r  can  be  expressed  as  functions  of 
the  element  node  displacements,  i.e.. 


uuk  -  Qp  u,'. 


(15) 


where  y=l...£2  and  |3=(£2+l),...,r. 

Substitution  of  (15)  into  (14)  thus  results  in  the  expression 


c l  a/  u;  Upk  =  Bf  u;  U  k  *  2  B?  Q|  Ut*  Qp"  Q*  U,1  Up' 


TP 


a 
k  b 


iir  nJP 


(16) 


where  y,p=l,...,£2  and  (3,co=(£2+l),...,r. 

If  we  now  take  the  derivative  of  (16)  with  respect  to  the  displacement  of  an  arbitrary  element  node 
a  in  degree-of-ffeedom  m,  we  obtain  3  x  £2  equations: 

Upk,  (17) 

where  a,p=l,...,£2  and  (3,co=(£2+l),...,r. 


rjl  Aap  uk  - 

cmk  Ajl  up  ” 


“H  n‘P  ^  BpP  <vta 


ml 


% 


kl 


Q£ 


BP®  r»j«  n1? 
jl 
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Since  the  Up  are  independent  of  each  other,  we  obtain  the  relations 


rjl  a*5  —  r^p  +  r  iP  nmp  +  rpP  n”r  +  rP®  nlY 

^ik  ^jl  ®ik  “im  ^(3k  ^(3i  °lm  ^(ok  * 


(18) 


where  y,p=l,...,£2  and  p,o=(£2+l),...,r. 


Recalling  that  from  the  transformation  (10)  we  have, 


ax1  asY  Yi 

i  ’ 


dq m  dqm 


(19) 


and  using  the  relations 


ax 1  aq  m  as Y  3q 


3q  m  dx ]  dqm  dx 


m  __i  -i 
jXy  ~  Oj, 


which  are  a  consequence  of  the  chain  rule  for  partial  differentiation,  we  can  multiply  expression  (13a) 
by  Xy1  Xpk  to  obtain 


A*5  xj  X  k  =  f-  —  iC  hi  xj  X  k  dv  =  V  8-  Sf, 
ji  *Y  P  Jv  dqm  aqn  aXJ  aXl  y  P  1 


(20) 


where  V  is  the  total  volume  of  the  element,  and  8j  is  the  Kronecker  delta. 


Multiplying  expression  (18)  by  xj  Xps  thus  results  in 


pjl  AYPYrYs-PISV  — 
'-'ik  Ajl  Ap  "  '■'ik  v  “ 


B?  *  b£  q£>  ♦  b£5  ♦  b£-  Qp?  Qj]  x;  x;,  (21) 
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where  y,p=l,...,£2  and  |3=(£2+l),...,r,  or,  upon  dividing  by  the  element  volume,  we  arrive  at  an  expression 
for  the  equivalent  element  material  stiffness  tensor;  i.e.. 


pji.  1  [  ryp  .  rtP  nmp  .  RpP  nmY  +  nmY  n"P  1  y i  y  1 

''ik  ~  "ik  +  “im  ^pk  +  “km  ^pi  +  “nm  ^pi  'vnk  * 

V 


(22) 


where  y,p=l,...,£2  and  |3=(£2+l),...,r. 

kv 

It  remains  to  develop  the  coefficients  Qpj  in  equation  (15),  which  expresses  the  displacements  of  the 
comer  points  (£2+1)  through  r  in  terms  of  the  displacements  of  the  element  nodes. 

For  the  A  comer  points  lying  on  the  surface  of  the  element,  we  require  that  the  net  force  on  these 
points  in  any  direction  tangent  to  the  surface  be  set  to  zero,  and  also  that  the  point  remains  on  the  surface 
as  the  element  deforms.  The  transformation  relations  between  global  cartesian  coordinates  (x1,  x2,  x3)  and 
isoparametric  coordinates  (q1,  q2,  q3)  must  be  developed  to  incorporate  these  constraints. 

At  the  location  of  one  of  the  comer  points  v  ((£2+1)  <  v  <  (£2+A)),  an  arbitrary  differential  length 
vector  can  be  expressed  in  the  global  cartesian  frame  as 

ds  =  dx  1  ij,  (23) 

where  Ij,  f2,  i3  correspond  to  the  units  vectors  along  the  global  cartesian  coordinate  directions.  From 
the  transformation  (10)  we  obtain 


dx 


(24) 


so  that 


ds  =  dq  k  Ij. 
3qk 


(25) 
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We  define  the  basis  vectors  for  the  isoparametric  coordinates  as 


— * 
b 


k  “ 


9x 
9q  k 


(26) 


to  obtain 


ds  =  bkdq  k. 


(27) 


in  the  isoparametric  system.  From  this  we  obtain  the  second  rank  metric  tensor. 


Sid  =  *>k  *  bj  = 


9x  1  9x  1 

dq*  V 


(28) 


In  general,  the  basis  vectors  of  the  isoparametric  system  are  not  orthogonal  so  that  the  off-diagonal 
components  of  the  metric  tensor  are  not  zero.  The  basis  vectors  are  tangent  to  the  isoparametric 
coordinate  curves  at  the  point  v.  Since  the  element  surface  on  which  v  lies  corresponds  to  a  surface  where 
one  of  the  coordinates,  say  qk,  is  a  constant,  the  basis  vectors,  tangent  to  the  other  two  coordinate  curves 
at  v,  are  tangent  to  the  element  surface  at  v,  (i.e.,  (b  and  (bnj^  (m  *  k  and  n  *  k)  lie  along  tangents 
to  the  element  surface  at  v,  with  the  surface  corresponding  to  qk  =  constant  in  the  isoparametric  system). 

From  (28)  we  can  develop  unit  tangent  vectors  along  the  surface,  given  by 


where  the  line  under  the  repeated  indices  indicates  that  summation  has  been  suspended. 

We  can  develop  the  reciprocal  base  vectors  for  the  isoparametric  system  from  the  vector  cross 
products 
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(30) 


b1 


_L  b2  x  b3,  b2  =  _L  b3  x  bv 

4i  Ji 


b3  =  _L  b,  x  b,, 


where  g  =  det  gM. 

These  relations  can  be  written  more  compactly  as 


(31) 


bi  = 


i  ^  (g® 

2^g 


W*. 


(32) 


where  e*1”1  is  the  familiar  permutation  symbol,  i.e., 

£kmn  =  1  if  k,  m,  n  is  an  even  pennutation  of  1,  2,  3; 
e10™1  =  -1  if  k,  m,  n  is  an  odd  permuatation  of  1,  2,  3;  and 
e^™  =  0  if  any  two  indices  are  the  same. 


The  reciprocal  base  vector 


corresponding  to  qk  =  constant,  i.e., 


J(v) 


at  v  is  therefore  in  the  direction  normal  to  the  element  surface 


if 


the  surface  corresponds  to  q1  =  constant,  [  b  lies  along  the  normal  and  [s2](V)  and  [§3]^) 
are  tangents  to  the  surface; 


if 


the  surface  corresponds  to  q2  =  constant, 
are  tangents  to  the  surface;  and 


b  2](v)  lies  along  the  normal  and  63](v)  and[ej] 


(v) 


J(v) 


if 


the  surface  corresponds  to  q3  =  constant,  [b3](V)  lies  along  the  normal  and  and[e2](v) 

are  tangents  to  the  surface. 
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We  now  define  the  contravarient  second  rank  tensor 


g 


kl 


This  now  allows  us  to  express  die  unit  normal  at  the  point  v  as 


1 


(33) 


(34) 


The  definitions  (29)  and  (34)  allow  us  to  apply  the  necessary  constraints  for  the  comer  points  (£2+1) 
through  (£2+A).  The  surface  comer  points  must  first,  however,  be  categorized  according  to  whether  a 
point  lies  along  one  of  the  edges  of  the  element,  or  whether  it  lies  on  an  element  face.  Let  the  points 
(£2+1)  through  (£2+E),  fKF.<A,  correspond  to  the  surface  comer  points  lying  on  the  edge  of  the  element, 
and  the  points  (£2+E+l)  through  (£2+A)  correspond  to  the  remaining  comer  points  that  lie  on  the  element 
faces. 

For  surface  comer  points  lying  along  an  element  edge  assume  that  the  comer  point  v  lies  on  an 
element  edge  corresponding  to  the  coordinate  curve  qr.  From  expression  (32),  the  reciprocal  base  vectors 
at  point  v  are  defined  as 


1  kmn  -  - 

-7=  e  (bm)(v)  x  (bn>(v)- 

2>/g 


(35) 


We  require  the  displacement  components  at  v  to  be  constrained  along  the  direction  of  bk(k*r)  such 
that  the  point  remains  on  the  surface  qk  =  constant,  corresponding  to  each  of  the  two  element  faces  that 
intersect  at  the  edge.  The  constraint  is  thus  expressed  as 

Uv*  Ii  *  (6  k)(v>  =  SP(q1,q2,q3)(v)  UpJ  •  (6  (36) 
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where  p=l,...,Q  and  v=(Q+l),...,(Q+E),  with  (e  k)(v)  representing  the  unit  vector  along  the  direction  of 
(bk)(v)  as  defined  in  (34). 

-  dx  * 

Since  bm  = _ L,  we  can  write 

dqm 


bk  =  _L  e’0™1  ill  —  t  x  L 
2sfe  aqm  aqn  J 


(37) 


or 


bk  -  _L  E1”  %  id  id  i„ 


2l/g 


5q  m  3q  n 


(38) 


where  we  have  used  the  expression  for  the  vector  cross  product 


h  x  Jj  =  eijl  II- 


(39) 


The  displacement  constraints  at  point  v  lying  on  an  element  edge  corresponding  to  the  isoparametric 
coordinate  curve  (f  can  thus  be  written  as 


J _ gknin  dx  1  dx3  TT  1 


2i/gVg“ 


kk  J  3q  m  dq  n 


u*  = 


1  gkmn  „  dx  1  dx 1 


2VgVg- 


«si  fir s  up ,  (40) 

kk  3q  m  5q  n 


or 


6*“™  —  —  Uvl  =  ekmn  ey!  ill  ill  SP(q1,q2,q3)(v)  U L1, 

3  aqm  3qn  3  3qm  9qn  P 


(41) 


where  k  *  r,  p=l,...,£l,  and  v=(£2+l ),..., (fl+E),  and  where  it  is  understood  that  the  expression  is  to  be 
evaluated  at  the  point  v. 
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The  third  constraint  expression  is  provided  by  the  requirement  for  zero  force  along  the  coordinate 
curve  qr  at  v.  To  impose  this  condition,  we  first  write  (8)  using  the  definition  (13b),  i.e.. 


FiV  -  BitV  *  B~U”  B^U*  *  B^Uck, 


(42) 


where  p=l,...,£2,  v,  (0=(f2+l),...,(ft+E),  p=(f2+E+l),...,(£2+A),  and  £=(£2+A+l),...,r. 


From  (26),  the  basis  vector  tangent  to  the  coordinate  curve  qr  is  expressed  as 


(43) 


The  condition  for  zero  force  along  q1  can  thus  be  enforced  by  forming  the  vector  dot  product  of 
(42)  with  (43)  and  setting  the  result  to  zero,  i.e.. 


F'V*‘  •  (B4>  ■  «• 


(44) 


or. 


3x  «  vco  t t  l  dx  p  vji  j t  l  3x  p  vC  tt  l  3x 

—  Ba  0.4  —  B,  U„  *^B«  U?=-^7 


RVP  TT  1 
Bil  Up » 


(45) 


r  corresponds  to  coordinate  curve  qr  along  the  element  edge  where  p=l,...,£2,  v,eo=(£2+l),...,(£2+E), 
p=(fl+E+l) . (Q+A),  and  C=(Q+A+1) . r. 

For  surface  comer  points  lying  on  an  element  face,  in  the  case  where  point  v  represents  a  comer  point 
lying  on  the  element  face  corresponding  to  cf  =  constant,  expression  (41)  provides  only  one  constraint 
equation,  restricting  the  displacement  components  of  point  v  such  that  it  remains  on  the  surface  whose 
normal  is  (6  k)(v)i.e„ 
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e*™  fiijj  ill  ili  Uv*  =  e*™  Ey!  ill  iiii  SPCqSq^q3)  U1, 

J  3q m  dq n  J  dqm  3q n  W  P 


(46) 


where  k  corresponds  to  qk  =  constant,  p=l,...,£2,  and  v=(Q+E+l),...,(£2+A). 

The  other  two  constraint  equations  are  obtained  from  the  requirement  for  zero  force  along  the 
directions  tangent  to  the  element  face  at  point  v,  i.e., 


^  -  °-  **  ^  •  MS  -  ®. 


(47) 


with 


•4S-HS 


defined  in  (29)  and  corresponding  to  unit  surface  tangent  vectors  at  point  v  on  the 


element  surface  representing  qk  =  constant 


These  conditions  lead  to  two  equations  of  the  form  (45), 


3x  1  „  v©  tt  1  .  1  R  VP  ti  1  .  1  R  v£  TI  i 

_  B„  U„  .  _  Bfl  U„  .  B, 


U?  *- 


3x  1  RvP  TI  l 

—  Up. 


(48) 


where  r*k,  p=l,...,Q,  co=(f2+l),...,(£i+E),  v,p=(Q+E+l),...,(f2+A),  ^=(fl+A+l),...,r,  and  where  it  is  again 
understood  that  (46)  and  (48)  are  to  be  evaluated  at  the  point  v. 

For  the  remaining  comer  points  ((Q+A+l)  through  T)  that  lie  internally  within  the  element’s 
boundaries,  we  require  zero  net  force  in  all  directions.  Expression  (42),  with  the  left-hand  side  set  to  zero, 
provides  the  constraint  equations  for  these  internal  points,  resulting  in 


B;ru,r  +  B;^u.r  +  =  -  b^u, 


V^TT  k  j.  R  V^TT  k 


W 


ik 


ik 


t  vPtt  k 


'ik  » 


(49) 


where  i=l,2,3,  (p=l,. 


„£2),  (co=(fl+l),...,(Q+E)),  (p=(Q+E+l),...,(Q+A)),  and  (v£=(Q+A+l) . T). 
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Expressions  (41),  (45),  (46),  (48),  and  (49)  together  form  a  system  of  3x(T-Q)  equations  in  the 
3x(T-£2)  unknown  displacements  for  the  A  surface  comer  points  and  the  (F-Q-A)  internal  comer  points. 
This  system  of  equations  can  be  solved  to  obtain  the  relations  (15),  expressing  the  displacement 
components  of  all  non-nodel  comer  points  in  terms  of  the  displacement  components  of  the  elements  nodes. 

Wedge  Prismatic  Elements 

The  constraint  equations  developed  above  are  suitable  for  rectangular  prismatic  elements  where  the 
edges  align  with  the  isoparametric  coordinate  curves.  For  wedge  prismatic  elements,  an  alternate  set  of 
constraint  equations  must  be  employed  for  comer  points  lying  on  the  element  surface. 

We  shall  assume  that  the  isoparametric  representation  of  the  wedge  prismatic  element  is  such  that  the 
two  triangular  faces  correspond  to  surfaces  of  q3  =  constant,  and  that  two  adjacent  four-noded  faces  of 
the  element  correspond  with  the  surfaces  q1  =  0  and  q2  =  0,  in  a  manner  that  forms  a  right-handed  frame 
of  reference.  The  intersection  of  the  third  four-noded  face  with  a  surface  represented  by  q3  =  constant 
results  in  a  curve  represented  by  q2  =  1  -  q1.  A  vector  lying  tangent  to  this  curve  at  a  point  v  can  be 
written  in  terms  of  the  basis  vectors  as 


(Q(V)  =  0>l)(v)  ~  tf^v)  = 


(  •  A 

dx  1  9x  1 


9q*  9q2 


(50) 


>) 


(C)(r)  and  the  basis  vector  (b3)(V)  thus  represent  two  independent  vectors  that  are  tangent  to  the  third 
four-noded  element  face  at  the  point  v. 

The  normal  to  this  element  face  at  point  v  can  be  obtained  from  the  cross  product 


N(v)  -  (C)(V)  *  (b^y)  = 


dx  1  _  ax1 
9q*  9q2 


Ii  x 


f 

axJ 


)(y) 


9q- 


(51) 


>) 


18 


or 


-  eijk 


{ dx  1 


f  \ 


ill 


(52) 


The  constraint  equations  for  wedge  prismatic  elements  can  now  be  developed.  As  with  the  rectangular 
prismatic  elements,  we  must  distinguish  between  surface  comer  points  that  lie  on  an  element  edge 
(i.e.,  points  (£2+1)  through  (£2+E»,  and  surface  comer  points  that  lie  on  element  faces  ((£2+E+l)  through 
(£2+A)). 

For  surface  comer  points  lying  along  an  element  edge,  the  two  displacement  constraint  equations 
provided  by  (41),  and  the  force  constraint  equation  provided  by  (45)  are  applicable  to  those  surface  comer 
points  lying  on  wedge  prismatic  element  edges  that  are  aligned  with  the  q1,  q2,  and  q3  coordinate  curves. 
The  only  two  element  edges  for  which  the  constraint  equations  (41)  and  (45)  do  not  apply,  are  the  edges 
corresponding  to  q2  =  1  -  q1  on  the  two  triangular  faces  of  the  element 

Since  the  triangular  faces  correspond  to  surfaces  q3  =  constant,  the  normal  to  these  faces  at  a  point 
v  is  obtained  from  the  reciprocal  base  vector  given  by  (38),  i.e.. 


<S3)(„  - 


2l/g~ 


e3mn  e- 


(  .\ 
dx  1  dx  1 


ijk 


aqm  aq"jv) 


(53) 


The  constraint  equations  for  a  point  v  lying  on  one  of  the  edges  corresponding  to  q2  =  1  -  q1  are  thus 
obtained  as  follows  along  the  direction  N^,  we  enforce  the  displacement  constraint 


•UI1i1=N2,-SP(q1,q2,q3)(l)Up1i1, 


(54) 
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where  p=l,...,£2  and  v=(£2+l),...,(£2+E),  so  that, 


A 


3x  1  _  dx  1 
3q*  3q2 


f  A 
3xJ 


;(v)^q  ^v) 


Uv  = 


f  .  A 

ax1  ax1 


aq1  3q2 


f  A 

3x J 


;<v)^  ;(v) 


sp(q1»q2.  q3)(v>  Up1,  (55) 


where  p=l,...,£2  and  v=(£2+l),...,(fl+E). 


Along  the  direction  (b3)(vj,  the  required  displacement  constraint  is 


(b3)(v)  •  Uv!  !j  =  (b3)v  •  S  p(q q2,  q3)(v)  Up1  llf 


(56) 


where  p=l,...,£2  and  v=(Q+l),...,(£2+E);  or. 


e3mn  ejjjl 


ax 1  axj 


aq  m  aq n 


Uv'  =  e3™  e.j, 


/V) 


3x 


A 


axJ 


dq  m  3q  n 


SP(q1,q2,q3)(l)  Up1,  (57) 


J(V) 


where  p=l,...,£2,  and  v=(£2+l),...,(£2+E),  and  along  the  direction  (C)(v),  which  is  tangent  to  the  edge  at 
v,  we  require  the  force  to  be  zero,  i.e.. 


V/v  _ 

•  F,-Ii  =  0, 


(58) 


resulting  in 


fax1 


L™  II  1  +  II  1 


BJCU 


dx 1  Rvp 
T — T  Bil 

^  >) 


(59) 


where  p=l,...,£2,  v,m=(Q+l),...,(f2+E),  ji=(£2+E+1),...,(£2+A),  and  ^=(i^+A+l),...,r. 
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For  surface  comer  points  lying  on  an  element  face,  only  those  comer  points  lying  on  the  element  face 
whose  normal  is  given  by  (52)  require  constraint  equations  that  are  different  from  (46)  and  (48). 


Expressions  (55)  and  (59)  provide  two  of  the  constraint  equations,  except  with  v  corresponding  to 
points  (Q+E+l)  through  (£2+A).  The  third  constraint  equation  corresponds  to  the  requirement  for  zero 

force  along  the  direction  of  the  basis  vector  (b3)(r)  that  can  be  obtained  from  (48)  by  setting  r  =  3,  i.e.. 


3x 1  D  vo  TT  l  9x 

9^  2  *+dq: 


b7  <  - 


9x  1  BvC  TT  l 


3q: 


UC  - 


3x 

*r 


Bavp  Up1, 


(60) 


where  p=l,...,£2,  co=(£2+l),...,(£2+E),  v,p=(f2+E+l),...,(fi+A),  and  £=(£2+A+l),...,r. 

Expressions  (55),  (57),  and  (59)  for  surface  comer  points  lying  on  the  edges  corresponding  to 
q2=l  -  q1,  and  (55),  (59),  and  (60)  for  surface  comer  points  lying  on  the  element  face  that  is  not  aligned 
with  the  isoparametric  coordinate  surfaces,  can  be  used  with  wedge  prismatic  elements  to  obtain  three 
equations  for  the  unknown  displacement  components  of  the  surface  comer  points.  These  relations, 
together  with  the  equations  provided  in  (41),  (45),  (46),  (48),  and  (49)  for  the  other  non-nodal  comer 
points,  can  be  grouped  to  form  the  3x(F-f2)  equations  required  for  obtaining  the  relations  (15),  when 
wedge  prismatic  elements  are  incorporated. 

3.  RESULTS 

In  this  section,  the  effective  stiffness  constants,  C^,  of  a  four-layered  cubic  element  with  rectangular 

faces  are  calculated  to  demonstrate  the  capability  of  this  newly  developed  formulation.  Results  are 
compared  to  those  calculated  by  Chou’s  model  showing  a  significant  difference  in  transverse  shear 
properties.  The  effects  of  stacking  sequence  of  layer  construction  on  the  effective  properties  of  die 
element  will  also  be  illustrated  and  discussed  in  detail. 

Figure  4  illustrates  the  coordinate  system  and  constitutive  relation  in  a  four-layer  laminated  block 
(0.2  x  0.2  x  0.2  in).  Coordinates  1  and  2  are  on  the  plane  of  laminate  plane.  The  ply  orientation  is 
defined  as  the  angle  between  fiber  direction  and  coordinate  1.  A  0°  ply  has  fibers  oriented  along 
coordinate  1.  The  effective  stiffness  components  are  illustrated  in  a  contracted  notation  (Cy,  i  and 
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Figure  4.  Constitutive  relation  for  an  anisotropic  four-layered  element. 


j  =  1  -  6)  for  convenience.  Two  layup  constructions,  a  cross-ply  [0/0/90/90]  and  an  angle-ply  [0/0/45/45], 
were  investigated.  Each  ply  has  an  equal  thickness  of  0.05  in.  In  fact,  each  ply  is  composed  of  10  unit 
directional  fiber  layers  with  thickness  of  0.005  in.  The  effective  properties  were  calculated  based  on 
IM7  graphite/8551  epoxy  material  whose  properties  are  shown  in  Table  1.  An  8-node  element  which 
utilizes  linear  transform  functions  in  Equation  4  for  analysis  was  used  to  calculate  effective  properties  of 
the  material  block. 
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Table  1.  Material  Properties  of  IM7/8551 


Ell 

= 

22.50E6  PSI 

E22 

= 

1.20E6  PSI 

E33 

= 

1.20E6  PSI 

vl2 

= 

0.33 

vl3 

= 

0.33 

v23 

= 

0.31 

G12 

= 

0.70E6  PSI 

G13 

- 

0.70E6  PSI 

G23 

= 

0.53E6  PSI 

3.1  Transverse  Shear  Properties.  Tables  2  and  3  show  the  comparison  of  effective  stiffness  for  both 
layup  constructions  [0/0/90/90]  and  0/0/45/45],  respectively.  Significant  differences  on  transverse  shear 
properties  and  shear  coupling  terms  were  found  for  both  cases.  The  transverse  shear  properties  (C^  and 
C55)  for  a  layup  construction  of  [0/0/90/90]  are  36%  lower  than  calculated  by  the  new  model.  The 
transverse  shear  properties  from  Chou’s  model  are  basically  calculated  from  a  volume  average.  His  model 
does  not  account  for  either  continuity  or  compatibility  of  the  materials  through  the  thickness.  A  linear 
deformation  was  also  made  by  Chou’s  model.  Therefore,  the  in-plane  properties  (Cjj  ,  =  1, 2,  and  6)  and 

transverse  normal  (C33)  were  found  to  be  identical.  These  properties  are  considered  to  be  exact  under  the 
assumption. 

For  an  angle-ply  layup  construction  [0/0/45/45],  transverse  shear  properties,  C44  and  C55,  are  different 
by  35%  and  40%,  respectively.  The  transverse  shear  coupling  terms  (C45  =  C54)  are  85%  of  difference. 
These  results  further  illustrate  the  importance  of  the  current  model.  In  fact,  larger  errors  may  be  obtained 
for  an  element  with  more  complex  ply  orientations,  stacking  sequence,  or  various  ply  thicknesses  by  using 
the  "volume  average"  approach. 

As  discussed  previously,  both  "volume  average"  and  "plate  theory"  approaches  cannot  accurately 
calculate  effective  properties  in  the  transverse  direction.  For  a  thick-section  structure,  the  transverse  shear 
properties  are  even  more  important  since  these  structures  generally  carry  more  shear  loads  than  thin-shelled 
structures.  For  finite  element  applications,  accurate  transverse  shear  properties  are  especially  important 
since  the  elements  are  3-D  blocks  with  arbitrary  shapes. 
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Table  2.  Comparison  of  the  Effective  Properties  of  Cross-Ply  Laminates 


Chou’s  Model 

0.1211E+08 

0.5883E+06 

0.505 1E+06 

0 

0 

0 

0.5883E+06 

0.1211E+08 

0.505 1E+06 

0 

0 

0 

0.505 1E+06 

0.505 1E+06 

0.1342E+07 

0 

0 

0 

0 

0 

0 

0.6033E+06 

0 

0 

0 

0 

0 

0 

0.6033E+06 

0 

0 

0 

0 

0 

0 

0.7000E+06 

1  3-D  Solid  Element  Model 

0.1211E+08 

0.5883E+06 

0.505 1E+06 

0 

0 

0 

0.5883E+06 

0.121 1E+08 

0.505 1E+06 

0 

0 

0 

0.5051E+06 

0.505 1E+06 

0.1342E+07 

0 

0 

0 

0 

0 

0 

0.4413E+06 

0 

0 

0 

0 

0 

0 

0.4413E+06 

0 

0 

0 

0 

0 

0 

0.7000E+06 

Table  3.  Comparison  of  the  Effective  Properties  of  Angle-Ply  Laminates 


Chou’s  Model 

0.1497E+08 

0.3117E+07 

0.5444E+06 

0 

0 

0.2694E+07 

0.3117E+07 

0.4194E+07 

0.465 8E+06 

0 

0 

0.2692E+07 

0.5444E+06 

0.465 8E+06 

0.1342E+07 

0 

0 

0.3933E+05 

0 

0 

0 

0.5670E+06 

0.4209E+05 

0 

0 

0 

0 

0.4209E+05 

0.6512E+06 

0 

0.2694E+07 

0.2692E+07 

0.3933E+05 

0 

0 

0.323 1E+07 

3-D  Solid  Element  Model 

0.1497E+08 

0.3117E+07 

0.5444E+06 

0 

0 

0.2694E+07 

0.3117E+07 

0.4194E+07 

0.4658E+06 

0 

0 

0.2692E+07 

0.5444E+06 

0.4658E+06 

0.1342E+07 

0 

0 

0.3933E+05 

0 

0 

0 

0.4182E+06 

0.227 1E+05 

0 

0 

0 

0 

0.227 1E+05 

0.4637E+06 

0 

0.2694E+07 

0.2692E+07 

0.3933E+05 

0 

0 

0.323 1E+07 
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3.2  Effects  of  Stacking  Sequence.  In  general,  the  plate  theory  assumes  constant  transverse  shear 
stress  distribution  through  the  thickness.  Accordingly,  the  effective  transverse  shear  constants  of  a 
laminate  calculated  from  the  plate  theory  approach  are  independent  of  the  stacking  sequence.  Recently, 
Roy  and  Kim  (1989)  showed  the  effects  of  stacking  sequence  on  transverse  shear  properties 
experimentally.  Models  based  on  the  deformations  of  a  beam  and  a  ring  subjected  to  specific  loading 
conditions  were  proposed  by  Roy  and  Tsai  (1992).  Their  model  reported  the  dependence  of  transverse 
shear  properties  on  stacking  sequence.  However,  only  two  specific  geometries  (beam  and  ring)  and 
loading  conditions  were  considered,  and  cannot  be  applied  to  a  generalized  case. 

Tables  4  and  5  illustrate  the  variations  of  electric  constants  in  cross-ply  [0/90]  and  angle-ply  [0/45] 
laminates  as  functions  of  stacking  sequence,  respectively.  In  the  cross-ply  laminate  case,  the  shear  elastic 
constant,  C44,  which  corresponds  to  shear  stress  and  strain  in  the  2-3  direction  (i.e.,  t23  and  y23)  increases 
as  the  90°  plies  are  located  away  from  the  laminate’s  midplane.  The  90°  plies  have  fibers  oriented  along 
coordinate  2  and  provide  more  shear  stiffness  in  the  2-3  direction.  Thus,  the  maximum  shear  stiffness, 
C44,  occurs  for  the  stacking  sequence  of  [90/0/0/90].  On  the  contrary,  the  shear  elastic  constant,  C55, 
which  corresponding  to  shear  stress  and  strain  in  the  1-3  direction  (xl3  and  yl3)  decreases  as  the  90°  plies 
are  moved  away  from  laminate’s  midplane.  Since  the  0°  layers  give  a  higher  shear  stiffness  in  this 
particular  direction,  the  laminate  with  layup  construction  of  [90/0/0/90]  has  the  smallest  shear  constant, 
C55. 

Similar  effects  of  stacking  sequence  on  transverse  shear  constants,  C44  and  C55»  and  shear  couplings, 
C45  and  C54,  were  observed  for  laminates  with  angle-ply  [0/45]  construction  shown  in  Table  5.  The  same 
conclusions  can  be  drawn  as  those  discussed  in  the  previous  section  for  the  cross-ply  laminates.  The  shear 
elastic  constant,  0*4'  which  corresponds  to  shear  stress  and  strain  in  the  2-3  direction  (i.e.,  x23  and  yl3), 
decreases  as  the  45°  plies  are  moved  away  from  the  laminate’s  midplane.  In  general,  the  effects  of 
stacking  sequence  are  more  significant  for  an  element  with  irregular  cross  sections,  complex  layup 
constructions,  and  nonsymmetric  stacking  sequence. 

4.  FINITE  ELEMENT  APPLICATIONS 

A  finite  element  preprocessor  was  developed  to  generate  finite  element  meshes  with  specific 
geometries.  The  effective  properties  of  each  element  are  calculated  individually  based  on  die  model. 
Figure  5  shows  the  finite  element  model  of  a  composite  rocket  motor  case  generated  by  the  preprocessor. 
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Table  4.  Effects  of  Stacking  Sequence  on  Transverse  Shear  Properties  of  Cross-Ply  Laminates 


[0/90/90/0] 

0.1211E+08 

0.5883E+06 

0.505 1E+06 

0 

0 

0 

0.5883E+06 

0.1211E+08 

0.505 1E+06 

0 

0 

0 

0.505 1E+06 

0.505 1E+06 

0.1342E+07 

0 

0 

0 

0 

0 

0  0.4278E+06 

0 

0 

0 

0 

0 

0 

0.4556E+06 

0 

0 

0 

0 

0 

0 

0.7000E+06 

[0/0/90/90] 

0.1211E+08 

0.5883E+06 

0.505 1E+06 

0 

0 

0 

0.5883E+06 

0.1211E+08 

0.505 1E+06 

0 

0 

0 

0.505 1E+06 

0.505 1E+06 

0.1342E+07 

0 

0 

0 

0 

0 

0  0.441 3E+06 

0 

0 

0 

0 

0 

0 

0.4413E+06 

0 

0 

0 

0 

0 

0 

0.7000E+06 

[90/0/0/90] 

0.121 1E+08 

0.5883E+06 

0.505 1E+06 

0 

0 

0 

0.5883E+06 

0.121 1E+08 

0.505 1E+06 

0 

0 

0 

0.505 1E+06 

0.505 1E+07 

0.1342E+07 

0 

0 

0 

0 

0 

0  0.4556E+06 

0 

0 

0 

0 

0 

0 

0.427 8E+06 

0 

0 

0 

0 

0 

0 

0.7000E+06 
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Table  5.  Effects  of  Stacking  Sequence  on  Transverse  Shear  Properties  of  Angle-Ply  Laminates 


[0/45/45/0] 

0.1497E+08 

0.3117E+07 

0.5444E+06 

0 

0 

0.2694E+07 

0.3117E+07 

0.4194E+07 

0.465 8E+06 

0 

0 

0.2692E+07 

0.5444E+06 

0.465 8E+06 

0.1342E+07 

0 

0 

0.3933E+05 

0 

0 

0 

0.4114E+06 

0.1576E+05 

0 

0 

0 

0 

0.3138E+05 

0.5925E+06 

0 

0.2694E+07 

0.2692E+07 

0.3993E+05 

0 

0 

0.3231E+07 

[0/0/45/45] 

0.1497E+08 

0.3117E+07 

0.5444E+06 

0 

0 

0.2694E+07 

0.3117E+07 

0.4194E+07 

0.4658E+06 

0 

0 

0.2692E+07 

0.5444E+06 

0.465 8E+06 

0.1342E+07 

0 

0 

0.3933E+05 

0 

0 

0 

0.5167E+06 

0.3441E+05 

0 

0 

0 

0 

0.3441E+05 

0.588 1E+06 

0 

0.2694E+07 

0.2692E+07 

0.3933E+05 

0 

0 

0.3231E+07 

[45/0/0/45] 

0.1497E+08 

0.3117E+07 

0.5444E+06 

0 

0 

0.2694E+07 

0.3117E+07 

0.4194E+07 

0.465 8E+06 

0 

0 

0.2692E+07 

0.5444E+06 

0.4658E+06 

0.1342E+07 

0 

0 

0.3933E+05 

0 

0 

0 

0.5204E+06 

0.3746E+05 

0 

0 

0 

0 

0.3746E+05 

0.5835E+06 

0 

0.2694E+07 

0.2692E+07 

0.3933E+05 

0 

0 

0.323 1E+07 
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Figure  5.  Finite  element  model  for  a  composite  rocket  motor  case. 

The  motor  case  consists  of  two  filament-wound  composite  components.  The  case’s  inner  region  is  a 
helically  wound  bottle  with  a  geodesic  winding  pattern.  The  outer  region  of  the  motor  case  is  a  cylinder 
with  a  cross-ply  layup  construction.  In  addition,  both  the  inner  and  outer  cases  were  constructed  with 
graphite  and  glass  composites. 


The  thickness  of  the  composite  varies  along  the  axial  direction  in  both  the  inner  and  outer  cases,  as 
shown  in  Figure  6.  Due  to  the  complexity  of  case  geometry,  the  general  elements  are  not  rectangular  and 
vaiy  along  the  arc  length.  The  elements  are  arbitrarily  shaped  and  contain  several  plies  with  various  fiber 
orientations  and  materials  through  the  thickness  in  some  areas.  In  fact,  it  is  typically  a  finite  element 
model  used  in  a  real-world  application;  the  capability  to  determine  the  effective  properties  in  the  proposed 
model  is  indeed  needed. 


For  a  geodesic  winding,  the  fibers  are  placed  along  the  shortest  path  on  the  case  surface  and  the 
winding  angle  (i.e.,  fiber  orientation)  varies  along  the  path  upon  the  geometry.  The  winding  angles  of  the 
innermost  layer  are  illustrated  and  plotted  along  the  arc  length  in  Figure  7.  The  winding  angles  vary 
dramatically  in  the  dome  and  nozzle  areas.  The  developed  preprocessor  calculates  the  winding  angle  at 
each  element  according  to  the  geodesic  path  and  case  geometry.  The  stiffness  components  (C^,)  of  the 

innermost  elements  along  the  arc  length  are  illustrated  in  Figure  8.  For  the  most  part,  effective  properties 
vary  significantly  in  most  materials.  The  variation  is  due  to  the  combined  effect  of  various  winding  angles 
and  layup  constructions. 

Clearly,  the  effective  properties  are  essential  to  achieving  an  accurate  finite  element  analysis. 
Currently,  all  the  commercial  packages  were  developed  using  either  "volume  average"  or  "laminated  plate 
theory"  to  determine  the  effective  properties.  In  general,  die  results  are  poor  for  structural  analyses  with 
complex  geometries  and  layup  constructions.  In  fact,  ABAQUS  suggests  that  no  skewed  elements  be  used 
in  the  finite  element  model  to  improve  the  accuracy  of  analyses. 
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Figure  8.  Stiffness  components  at  the  innermost  layer  vary  along  the  arc-length. 
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5.  CONCLUSIONS 


Based  on  strain  energy  approaches  and  finite  element  techniques,  an  effective  property  model  was 
developed  to  determine  the  properties  of  an  arbitrarily  shaped  element  with  multi-material  regions.  The 
model  is  especially  suitable  for  3-D  finite  element  application  due  to  the  accurate  transverse  shear 
properties.  The  effective  material  stiffness  calculated  by  the  model  was  compared  to  these  by  Chou’s 
model.  Significant  differences  in  transverse  shear  properties  were  found  for  a  four-ply  cubic  element. 
The  comparison  illustrated  the  lack  of  accuracy  of  currently  available  models.  Effects  of  stacking 
sequence  on  transverse  properties  were  identified  and  discussed  in  detail. 

Having  accurate  transverse  shear  properties  and  the  capability  to  model  arbitrarily  shaped  elements 
are  particularly  important  for  finite  element  applications,  especially  for  thick-section  composites  with  near- 
net  shape  geometries  subjected  to  complex  loadings.  A  preprocessor  was  developed  using  the  effective 
property  model  to  generate  properties  for  the  finite  element  model.  The  preprocessor  currently  has 
capabilities  to  generate  material  properties  for  a  finite  element  model  for  an  axisymmetric  filament-wound 
case  and  several  other  geometries  of  interest.  The  preprocessor  is  developed  to  be  used  with  DYNA3D 
and  ABAQUS  finite  element  codes  to  perform  dynamic  analysis  of  composite  structures  with  complex 
thick-section  geometry. 
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PO  BOX  12211 

RESEARCH  TRIANGLE  PARK  NC  27709-2211 

3  US  ARMY  RESEARCH  OFFICE 
ENGINEERING  SCIENCES  DIV 
ATTN  R  SINGLETON 

G  ANDERSON 
K  IYER 
PO  BOX  12211 

RESEARCH  TRIANGLE  PARK  NC  27709-2211 

2  PROJECT  MANAGER 
SADARM 

PICATINNY  ARSENAL  NJ  07806-5000 

2  PROJECT  MANAGER 

TANK  MAIN  ARMAMENT  SYSTEMS 
ATTN  SFAE  AR  TMA  COL  BREGARD 
KKIMKER 

PICATINNY  ARSENAL  NJ  07806-5000 

1  PROJECT  MANAGER 

TANK  MAIN  ARMAMENT  SYSTEMS 
ATTN  SFAE  AR  TMA  MD  R  KOWALSKI 
PICATINNY  ARSENAL  NJ  07806-5000 

2  PEO  FIELD  ARTILLERY  SYSTEMS 
ATTN  SFAE  FAS  PM 

D  ADAMS 
T  MCWILLIAMS 

PICATINNY  ARSENAL  NJ  07806-5000 

1  PEO  FIELD  ARTILLERY  SYSTEMS 
ATTN  SFAE  FAS  PM  H  GOLDMAN 
PICATINNY  ARSENAL  NJ  07806 

2  PROJECT  MANAGER  AFAS 
ATTN  G  DELCOCO 

J  SHIELDS 

PICATINNY  ARSENAL  NJ  07806-5000 
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2  NASA  LANGLEY  RESEARCH  CTR 
MS  266 

ATTN  AMSRL  VS 
WELBER 
F  BARTLETT  JR 
HAMPTON  VA  23681-0001 

2  COMMANDER 
DARPA 

ATTN  J  KELLY 
B  WILCOX 
3701  N  FAIRFAX  DR 
ARLINGTON  VA  22203-1714 

2  COMMANDER 

WRIGHT  PATTERSON  AIR  FORCE  BASE 
ATTN  WL  FTV  A  MAYER 
DAYTON  OH  45433 

2  NAVAL  SURFACE  WARFARE  CTR 
DAHLGREN  DIV 
CODE  G33 

DAHLGREN  VA  224488 

1  NAVAL  RESEARCH  LABORATORY 
CODE  6383 
ATTN  I WOLOCK 
WASHINGTON  DC  20375-5000 

1  OFFICE  OF  NAVAL  RESEARCH 

MECH  DIV  CODE  1132SM 
ATTN  YAPA  RAJAPAKSE 
ARLINGTON  VA  22217 

1  NAVAL  ORDNANCE  STATION 

ADVANCED  SYSTEMS  TECHNOLOGY  BR 
ATTN  D  HOLMES 
CODE  2011 

LOUISVILLE  KY  40214-5245 

1  DAVID  TAYLOR  RESEARCH  CTR 

SHIP  STRUCTURES  AND  PROTECTION  DEPT 
ATTN  J  CORRADO  CODE  1702 
BETHESDA  MD  20084 

2  DAVID  TAYLOR  RESEARCH  CTR 
ATTN  R  ROCKWELL 

W  PHYILLAIER 
BETHESDA  MD  20054-5000 


1  DEFENSE  NUCLEAR  AGENCY 

INNOVATIVE  CONCEPTS  DIV 
ATTN  DR  R  ROHR 
6801  TELEGRAPH  RD 
ALEXANDRIA  VA  22310-3398 

1  DR  FRANK  SHOUP 

EXPEDITIONARY  WARFARE  DIV  N85 
2000  NAVY  PENTAGON 
WASHINGTON  DC  20350-2000 

1  OFFICE  OF  NAVAL  RESEARCH 

ATTN  MR  DAVID  SIEGEL  351 
800  N  QUINCY  ST 
ARLINGTON  VA  22217-5660 

1  NAVAL  SURFACE  WARFARE  CTR 

ATTN  CODE  G30 
JOSEPH  H  FRANCIS 
DAHLGREN  VA  22448 

3  NAVAL  SURFACE  WARFARE  CTR 
ATTN  CODE  G33 

JOHN  FRAYSSE 
ROD  HUBBARD 
GIL  GRAFF 
DAHLGREN  VA  22448 

1  DEFENSE  NUCLEAR  AGENCY 

INNOVATIVE  CONCEPTS  DIV 
ATTN  LTC  JYUJI D  HEWITT 
6801  TELEGRAPH  RD 
ALEXANDRIA  VA  22448 

1  CDR  NAVAL  SEA  SYSTEMS  CMD 

ATTN  D  LIESE 

2531  JEFFERSON  DAVIS  HIGHWAY 
ARLINGTON  VA  22242-5160 

1  NAVAL  SURFACE  WARFARE 

ATTN  CODE  D4  MARY  E  LACY 
17320  DAHLGREN  RD 
DAHLGREN  VA  22448 

4  DIRLLNL 

ATTN  R  CHRISTENSEN 
S  DETERESA 
F  MAGNESS 
M  FINGER 
PO  BOX  808 
LIVERMORE  CA  94550 
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1  DIRLANL 

ATTN  D  RABERN 
MEE  13  MS  J  576 
PO  BOX  1633 
LOS  ALAMOS  NM  87545 

1  DIRLANL 

ATTN  J  REPPA  MS  F668 

PO  BOX  1663 

LOS  ALAMOS  NM  87545 

1  OAK  RIDGE  NATIONAL  LABORATORY 

ATTN  R  M  DAVIS 
PO  BOX  2008 

OAK  RIDGE  TN  37831-6195 

4  DIR  SANDIA  NATIONAL  LABORATORIES 

APPLIED  MECHANICS  DEPARTMENT 
ATTN  DIVISION  8241 
G  BENEDETTI 
KPERANO 
D  DAWSON 
PNIELAN 
PO  BOX  969 

LIVERMORE  CA  94550-0096 

1  BATTELLE 

ATTN  C  R  HARGREAVES 
505  KING  AVE 
COLUMBUS  OH  43201-2681 

1  PACIFIC  NORTHWEST  LABORATORY 

ATTN  M  SMITH 
PO  BOX  999 
RICHLAND  WA  99352 

1  DIRLLNL 

ATTN  M  MURPHY 
PO  BOX  808  L  282 
LIVERMORE  CA  94550 

1  PURDUE  UNIVERSITY 

SCHOOL  OF  AERO  AND  ASTRO 
ATTN  C  T  SUN 
W  LAFAYETTE  IN  47907-1282 

1  STANFORD  UNIVERSITY 

DEPARTMENT  OF  AERONAUTICS  AND 
AEROBALLISTICS  DURANT  BUILDING 
ATTN  S  TSAI 
STANFORD  CA  94305 


1  UCLA 

MANE  DEPT  ENGR  IV 
ATTN  H  THOMAS  HAHN 
LOS  ANGELES  CA  90024-1597 

2  U  OF  DAYTON  RESEARCH  INSTITUTE 
ATTN  RAN  Y  KIM 

AJTTKROY 

300  COLLEGE  PARK  AVE 
DAYTON  OH  45469-0168 

2  UNIVERSITY  OF  DELAWARE 

CTR  FOR  COMPOSITE  MATERIALS 
ATTN  J  GILLESPIE 
201  SPENCER  LABORATORY 
NEWARK  DE  19716 

1  THE  UNIVERSITY  OF  TEXAS  AT  AUSTIN 
CENTER  FOR  ELECTROMECHANICS 
ATTN  J  PRICE 

10100  BURNET  RD 
AUSTIN  TX  78758-4497 

2  VIRGINIA  POLYTECHNICAL  INSTITUTE 
AND  STATE  UNIVERSITY 

DEPT  OF  ESM 
ATTN  MICHAEL  W  HYER 
KENNETH  L  REEFSNIDER 
BLACKSBURG  VA  24061-0219 

1  AAI  CORPORATION 

PO  BOX  126 

HUNT  VALLEY  MD  21030-0126 

1  SAIC 

ATTN  DAN  DAKIN 
2200  POWELL  ST  STE  1090 
EMERYVILLE  CA  94608 

1  SAIC 

ATTN  MILES  PALMER 
2109  AIR  PARK  RD  S  E 
ALBUQUERQUE  NM  87106 

1  SAIC 

ATTN  ROBERT  ACEBAL 

1225  JOHNSON  FERRY  RD  STE  100 

MARIETTA  GA  30068 
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1  SAIC 

ATTN  DR  GEORGE  CHRYSSOMALLIS 
3800  W  80TH  STREET 
STE  1090 

BLOOMINGTON  MN  55431 


OUN  CORPORATION 
FLINCHBAUGH  DIV 
ATTN  E  STEINER 
B  STEWART 
PO  BOX  127 
RED  LION  PA  17356 


4  ALLIANT  TECHS YSTEMS  INC. 

ATTN  C  CANDLAND 
J  BODE 
R  BECKER 
K  WARD 
600  2ND  ST  NE 
HOPKINS  MN  55343-8367 


1  OUN  CORPORATION 
ATTN  L  WHITMORE 
10101  9TH  ST  NORTH 
ST  PETERSBURG  FL  33702 


1  CUSTOM  ANALYTICAL  ENGINEERING 
SYSTEMS  INC 
ATTN  A  ALEXANDER 
STAR  ROUTE  BOX  4A 
FLINTSTONE  MD  21530 


1  PROJECTILE  TECHNOLOGY  INC 

515  GILES  ST 

HAVRE  DE  GRACE  MD  21078 

1  ALLEN  BOUTZ 

NOESIS  INC 
1110  N  GLEBE  RD 
STE  250 

ARLINGTON  VA  22201-4795 

1  GENERAL  DYNAMICS 

LAND  SYSTEMS  DIVISION 
ATTN  D  BARTLE 
PO  BOX  1901 
WARREN  MI  48090 


4  INSTITUTE  FOR  ADVANCED  TECHNOLOGY 

ATTN  T  KIEHNE 
H  FAIR 
P  SULLIVAN 
IMCNAB 

4030  2  W  BRAKER  LN 
AUSTIN  TX  78759 

1  INTERFEROMETRICS  INC 

ATTN  R  LARRIVA  VICE  PRESIDENT 
8150  LEESBURG  PIKE 
VIENNA  VA  22100 
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ABERDEEN  PROVING  GROUND 

DIRUSARL 
ATTN  AMSRLSC 

W  MERMAGEN 
CWSTUREK 
AMSRL  SC  S  A  MARK 
AMSRL  MA  P  L  JOHNSON 
AMSRL  MA  PD  T  CHOU 
AMSRL  MA  PA  D  GRANVILLE 
AMSRL  MA  MA  G  HAGNAUER 
AMSRL  SL  B  P  DIETZ  328 
AMSRL  SL  BA  J  WALBERT  1065 
AMSRL  SL  BL  D  BELY  328 
AMSRL  SL  I D  HASKILL  1065 
AMSRL  WTP 
A  HORST  390A 
E  SCHMIDT  390A 
AMSRL  WT  PA 
T  MINOR  390 
D  KOOKER  390A 
AMSRL  WTPB 
P  PLOSTINS  390 
D  LYON 

AMSRL  WT  PC  R  FIFER  390A 
AMSRL  WT  PD 
B  BURNS  390 
W  DRYSDALE  390 
J  BENDER  390 
T  ERLINE  390 
D  HOPKINS  390 
S  WILKERSON  390 
R  KASTE  390 
L  BURTON  390 
J  TZENG  390 
R  LIEB  390 
G  GAZONAS  390 
C  HOPPEL  390 
AMSRL  WT  PD  ALC 
A  ABRAHAMIAN 
•  M  BERMAN 

A  FRYDMAN 
T  LI 

W  MCINTOSH 


ABERDEEN  PROVING  GROUND  (CONTINUED) 

AMSRL  WT  T  W  MORRISON  309 
AMSRL  WT  TA 
W  GILLICH  390 
W  BRUCHEY  390 
AMSRL  WT  TC 
K  KIMSEY  309 
R  COATES  309 
W  DE  ROSSET  309 
AMSRL  WTTD 
D  DIETRICH  309 
G  RANDERS  PEHRSON  309 
A  DAS  GUPTA  309 
J  SANTIAGO  309 
AMSRL  WT  W  C  MURPHY  120 
AMSRL  WT  WA 
H  ROGERS  394 
B  MOORE  394 
AMSRL  WT  WB 
F  BRANDON  120 
W  D  AMICO  120 

AMSRL  WT  WC  J  BORNSTEIN  120 
AMSRL  WTWD 
A  ZIELINSKI  120 
J  POWELL  120 
AMSRL  WT  WE 
J  LACETERA  120 
J  THOMAS  394 
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USER  EVALUATION  SHEET/CHANGE  OF  ADDRESS 

This  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the  reports  it  publishes.  Your  comments/answers 
to  the  items/questions  below  will  aid  us  in  our  efforts. 

1.  ARL  Report  Number  ARL-TR-1051 _ Date  of  Report  April  1996 _ 

2.  Date  Report  Received _ 

3.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related  project,  or  other  area  of  interest  for  which  the  report 

will  be  used.) _ 


4.  Specifically,  how  is  the  report  being  used?  (Information  source,  design  data,  procedure,  source  of  ideas,  etc.) 


5.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far  as  man-hours  or  dollars  saved,  operating  costs 
avoided,  or  efficiencies  achieved,  etc?  If  so,  please  elaborate. _ 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future  reports?  (Indicate  changes  to 
organization,  technical  content,  format,  etc.) _ 


Organization 


CURRENT  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 


City,  State,  Zip  Code 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the  Current  or  Correct  address  above  and  the 
Old  or  Incorrect  address  below. 


Organization 


OLD  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 


City,  State,  Zip  Code 


(Remove  this  sheet,  fold  as  indicated,  tape  closed,  and  mail.) 
(DO  NOT  STAPLE) 


DEPARTMENT  OF  THE  ARMY 


OFFICIAL  BUSINESS 


BUSINESS  REPLY  MAIL 

FIRST  CLASS  PERMIT  NO  0001  ,APG,MD 


POSTAGE  WILL  BE  PAID  BY  ADDRESSEE 


DIRECTOR 

U.S.  ARMY  RESEARCH  LABORATORY 
ATTN:  AMSRL-WT-PD 

ABERDEEN  PROVING  GROUND,  MD  21005-5066 
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